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Abstract 



We present a new complete set of Lagrangian relativistic hydrodynamical equations 
describing the transfer of energy and momentum between a standard fluid and a 
radiation fluid in a general non-stationary spherical flow. The new set of equations 
has been derived for a particular application to the study of the cosmological 
Quark-Hadron transition but can also be used in other contexts. 
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1. Introduction 

An important problem in the study of relativistic plasmas is the correct 
treatment of the coupling between radiation and matter. The assumptions of the 
classical theory of radiative transfer are no longer adequate near to astrophysical 
compact objects or for some cosmological situations and a relativistic theory of 
radiation hydrodynamics is required for these. 

In this paper we derive a new system of Lagrangian relativistic hydrodynamic 
equations for describing general non-stationary spherical flows in which a transfer 
of energy and momentum occurs between a standard perfect fluid and a generalized 
"radiation fluid" (i.e. a fluid composed of any effectively zero rest-mass particles 
having longer mean- free-path than those of the standard fluid). Our primary 
interest is in application of these equations to cosmological phase transitions and 
to spherical accretion onto black holes. This paper presents the basic system of 
equations; further details of the applications will be given in subsequent papers. 

Consistent treatments of relativistic radiation hydrodynamics have only been 
developed fairly recently and among these, particular interest has been focussed on 
the PSTF tensor formalism devised by K.S. Thorne in 1981 [1, 2] and the covariant 
flux-limited diffusion theory (FLDT) approach proposed by A.M. Anile et al. [3]. 
Both methods represent approximations to a problem which is too complicated to 
handle in full detail at the present time and there has been some debate about their 
respective merits. Since the PSTF method has already been used for a number of 
successful calculations [4, 5, 6], we have followed this approach here but we await 
with interest results from applications of the alternative method. We will comment 
in subsequent sections on our reasons for thinking that the approach which we are 
following is a satisfactory one in the present context. 

In a recent paper, M.-G. Park [7] has presented a system of relativistic ra- 
diation hydrodynamic equations for treating a time-dependent spherical accretion 
flow onto a black hole or neutron star. The approach there is rather different 
from that of the present work and, in particular, Park writes the fluid equations 
in an Eulerian (fixed) frame. While this form is convenient for some purposes, the 
Lagrangian formulation has great advantages for many practical applications in 
computations of one-dimensional time-dependent flows (having planar, cylindrical 
or spherical symmetry). 

In Section 2 we briefly describe the fundamentals of the PSTF formalism 
and in Sections 3 and 4 we discuss the derivation of the hydrodynamical equations 
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for the radiation fluid and the standard fluid. We here follow the approach of the 
earlier papers [8, 9, 10] which are relevant background for the present work. The 
subsequent sections introduce further additional equations which are necessary 
for treating the problem of bubble growth during the cosmological Quark-Hadron 
transition. Section 5 contains a short description of the scenario for the transition 
and in Section 6 we present interface junction conditions and characteristic equa- 
tions for the treatment of a discontinuity surface between two phases of the stan- 
dard fluid. Discussion of the numerical implementation of the formalism presented 
here will be delayed until a subsequent paper [11] but, as evidence that a successful 
implementation can be made, some sample results are shown in Section 7. 

In view of our application to a problem involving input from particle physics, 
we use in this paper a system of units in which c = h = k B = 1 rather than 
c = G = 1 as is usual for calculations in general relativity. This also has the 
advantage that the gravitational source terms are clearly identified because of 
retention of the constant G in the equations. (Note that Thorne [1, 2] adopted the 
slightly different convention c = h = 1). We use a space-like signature (— , +, +, +). 
Greek indices are taken to run from to 3 and Latin indices from 1 to 3. Covariant 
derivatives are denoted with a semi-colon and partial derivatives with a comma. 



2. The PSTF Tensor Formalism 

The PSTF tensor formalism is a technique for solving the general relativistic 
form of the radiative transfer equation, which describes the variation of the radi- 
ation field as it propagates through a standard fluid. The relativistic form of this 
is straightforward and can be written as 



where M is the distribution function for the photons (a relativistic invariant), / 
is a non-affine parameter measuring the proper spatial distance travelled by the 
photons as seen from the standard fluid, and E is a source function. Note that the 
total derivative is taken not just in the space-time but rather in the phase-space 
since M = M{x a ,p a ), where p is the photon four-momentum. 

The fundamental idea of the PSTF method consists of replacing equation 
(1) (which is in a concise form but embodies enormous complexity) by a hierarchy 
of moment equations written in terms of Projected Symmetric and Trace-Free 
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(PSTF) tensors which are suitably defined at each point in the projected tangent 
space to the fluid four- velocity u. The k-th moments of M and E are 



where v is the specific frequency under consideration, dV p is the invariant 
momentum-space volume element on the light cone and S(y) is the Dirac delta 
function. In the following, we will refer to JAJ* lm " ak and S l/ ai " ak simply as the 
k-th moment and source moment respectively. 

The expressions (2) and (3) can be integrated over frequency and a clear 
physical interpretation can be given for the first three integrated moments of the 
hierarchy: M. (the zero-th moment) is the energy density of the radiation, M a 
(the first moment) is the radiative energy flux, and M a/3 (the second moment) is 
the shear stress tensor of the radiation fluid (each quantity being measured in the 
local rest frame of the standard fluid). The stress-energy tensor for the radiation 

j g completely defined in terms of the first three moments and higher order 
moments do not enter into this definition. The expression for it is 



where P a/3 (= g a/3 + u a u^) is the projection operator orthogonal to u. A conse- 
quence of this is that, if the hierarchy is truncated at the second order, it is possible 
to derive the equations governing the hydrodynamics of the radiation fluid in a 
particularly simple way by starting from the conservation laws of energy and mo- 
mentum. If, on the other hand, orders higher than the second are retained, it 
is necessary to make direct use of the appropriate hierarchy of equations derived 
from (1). 

In the case of planar or spherical symmetry, the 2k + 1 independent compo- 
nents of each k rank tensor depend on a single scalar variable so that the tensor 
formalism reduces to a purely scalar one. This simplification has made it possible 
for the method to be used for a number of astrophysical applications. The hierar- 
chy of integrated scalar moment equations, into which equation (1) is recast, has 
the property that, for any k, the first k equations involve the first k + 1 moments. 
In order to use this scheme for making calculations, it is necessary to truncate 






(3) 




(4) 
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the moment hierarchy at some finite order by introducing a closure relation which 
specifies the value of the highest moment used in terms of lower ones and which 
is derived on the basis of physical considerations. 

In the next section we use the PSTF approach to derive a new set of hydro- 
dynamical equations describing the coupling between radiation and matter in the 
case of a non-stationary spherically symmetric flow, with the moment equation 
hierarchy being truncated at the second order. 

3. Hydrodynamics of the Radiation Fluid 

We consider a spherically symmetric flow in which a transfer of energy and 
momentum occurs between a generalized "radiation fluid" and a standard fluid. 
It is convenient to use a Lagrangian frame comoving with the standard fluid and 
the spherically symmetric line element 

ds 2 = -a 2 dt 2 + b 2 d/j, 2 + R 2 (d9 2 + sin 2 9 dip 2 ), (5) 

where \x is a comoving radial coordinate and R is an associated Eulerean coordinate 
(the Schwarzschild circumference coordinate). 

We here describe the radiation hydrodynamics using the first two moment 
equations together with a closure relation and the calculations then involve the first 
three moments, together with the first two source moments. The maximum errors 
in the calculated values of the radiation variables which arise when truncating 
at the second order, are typically ~ 15% [8]. In view of the fact that other 
uncertainties in the specification of the problem are of a comparable order, we 
regard this level of approximation as acceptable for the present purposes. Similar 
degrees of accuracy are reported for calculations using FLDT schemes [12]. 

For spherical symmetry the first three moments can be written as 

M = w , (6) 

M a = Wl e?, (7) 

and 

M ap = w 2 (e?4 - ±e£e£ - \e%ety , (8) 

where w$, w\ and w 2 are the scalar moments and (e 5 , e-, e^, e^) is the orthonor- 
mal tetrad carried by an observer comoving with the standard fluid. 
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The quantities wq, u>\ and 102 all have direct physical interpretations corre- 
sponding to those of the related tensor moments: wq and w\ are the energy density 
and flux of the radiation in the fluid rest frame, while W2 is the shear stress scalar 
of the radiation. The scalar source moments sq and s\, defined in a similar way to 
wq and wi, also have direct physical interpretations, representing the transfer of 
energy and momentum between the two fluids. We use the following expressions 
for these scalar source moments 

so = j(e-w ) + (s ) sc , (9) 

.. = -!£. do) 

where A is the effective mean-free-path of the radiation particles as they move 
through the standard fluid, (sq) sc is a term expressing the contribution due to 
scatterings, whose form depends on the specific problem, and e is the energy 
density for radiation in thermal equilibrium with the standard fluid. Assuming 
that it follows a black-body law, e can be written as 

< = (») 

with g R being the number of degrees of freedom of the radiation fluid and T F the 
temperature of the standard fluid. 

Having restricted ourselves to the use of the first three scalar moments and 
the first two scalar source moments, it is convenient to derive the radiation hydro- 
dynamical equations by means of the standard conservation laws for the energy 
and momentum of the radiation fluid applied to the stress-energy tensor (4). Fol- 
lowing this procedure, which is equivalent to using the scalar moment equations, 
we then write the following three radiation hydro dynamical equations 



-u a Tf.p = so, (12) 
njKTf^ = * (13) 
w 2 = f E w , (14) 



where n is a radial spacelike unit vector normal to u. 

The term f E in the closure relation (14) is a variable Eddington factor and 
indicates the degree of anisotropy of the radiation. It can take values ranging 
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from for complete isotropy (which could, for example, be caused by the medium 
being extremely optically thick) to 2/3 for complete anisotropy. A key point in 
the present technique is that an expression for f E has to be supplied, constructed 
on the basis of physical considerations and how this is done is, to some extent, 
ad hoc. However, experience has shown that as long as the expression has the 
correct asymptotic behaviour in any relevant limits, results do not normally de- 
pend sensitively on the precise form used as long as it gives a suitably smooth join 
between the limits [4]. This is something which needs to be checked in any par- 
ticular application but provided that the outcome of such a check is satisfactory, 
it is reasonable to proceed with confidence. 

Making use of the stress-energy tensor (4) and of the line element (5), equa- 
tions (12) and (13) can be written explicitly as 

(15) 

, x . afl \ 4a iM , fb, t R,t\ a ( a M 3.R M \ 

M^ + l{-^o + w 2 j^^ Wo + 2^ T + -j Wl + -^- + ^j W2 = as 1 , 

(16) 

(see Appendix A for details). Equations (14) - (16) are our final form of the hy- 
drodynamical equations for the radiation fluid and they need to be solved together 
with the corresponding hydro dynamical equations for the combined fluids which 
will be discussed in the next section. 

4. Hydrodynamics of the Combined Fluids 

The derivation of the hydro dynamical equations for the combined fluids (i.e. 
standard fluid and radiation fluid) is, in principle, more straightforward. We 
follow the conventions and notation of the earlier papers [9, 10] (except that here 
we always write partial derivatives using a comma) and start from the Einstein 
field equations 

Ra/3 - -^9apR = 8nGT a p, (17) 

with T a/3 = T^ 3 + T£@ being the total stress-energy tensor, and treat the standard 
fluid as perfect so that T"^ = (e + p)u a u^ + pg af3 . The four independent equa- 
tions which follow from (17) can be used to write out explicitly the conservation 
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equations for energy and momentum of the combined fluids and the continuity 
equation for the standard fluid 



-u a T a % = 0, (18) 

njP'aT"^ = 0, (19) 
{pu a ), a = 0, (20) 

where p is the relative compression factor expressing the variation in the proper 
volume of co-moving elements of the standard fluid. (For a classical standard fluid 
composed of non-relativistic particles, p can represent the rest- mass density.) 

After some algebra, whose relevant steps are shown in Appendix A, it is 
possible to write the following set of hydrodynamical equations 

GM~ 



u* = -a 



R 2 



(21) 



R,t = au, (22) 
( P R%_ Ju,-4,GoR Wl \ 



pR 2 V R» 

e,t = wp, t - as , (24) 

Kg _ wp^-e^ + b Sl ^ 
aw pw 

ii \ 

(26) 



M,„ = 4tt# 2 # M + w + ^w^j 



2GM\ 1/2 1 



) =b R - (27) 
i= sW <28) 

Here u is the radial component of fluid four velocity in the associated Schwarzschild 
(Eulerian) frame, V is the general relativistic analogue of the Lorentz factor, and 
w is the specific enthalpy (w = (e + p) / p) . The generalized mass function M can 
also be calculated using the alternative equation 

M ;t = -4irR 2 R, t (p + ^w + + w 2 ^j . (29) 

The set of hydrodynamical equations (14) - (16) and (21) - (28) needs to 
be supplemented by an equation of state, relating the energy density e, the pres- 
sure p and temperature T of the standard fluid. The basic set of equations is 
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then completed and can be used for describing the transfer of energy and momen- 
tum between a standard fluid and a radiation fluid in a general non-stationary 
relativistic flow. 

It is worthy of note that the gravitational source terms appearing in the 
combined set of equations introduce only a very minor additional complication; 
for spherical symmetry, a general relativistic calculation is only marginally more 
complicated than a special relativistic one. 

We next turn to the introduction of the further equations necessary for treat- 
ing the problem of bubble growth during the cosmological Quark-Hadron tran- 
sition. Here, there are two different phases of the standard fluid separated by a 
phase interface (the bubble surface) and so the above system of equations (which 
continue to hold for the bulk of each phase) needs to be supplemented by a treat- 
ment of the interface. Before discussing how this is done (in Section 6), Section 5 
gives a brief introduction to the scenario within which we are working. 

5. The Cosmological Quark— Hadron Transition 

According to the standard Hot Big-Bang model, the early universe experi- 
enced a succession of phase transitions and breakings of symmetry which would 
have influenced the subsequent evolution. The last transition in this succes- 
sion (which includes Inflation and the Electro- Weak transition) is usually thought 
to have been the Quark-Hadron transition at which strongly interacting matter 
passed from being a plasma of unconfined quarks and gluons to a plasma in which 
the quarks and gluons were confined within hadrons. Since it is thought to have 
been the last of the cosmological transitions, any remnants which it left behind 
could have had particular significance and so it is of considerable interest to carry 
out a detailed hydrodynamical study in order to see how it would have proceeded. 

Chronologically, the transition is estimated to have taken place about 10/US 
after the Big-Bang, when the temperature of the universe was 200 MeV, the 
mean density was w 10 15 g cm -3 and the horizon scale was « 10 km. At present, 
QCD lattice calculations carried out to investigate the nature of the transition 
seem to indicate that it may well be a continuous one [13]. However, these calcu- 
lations are extremely complicated and the results are not yet definitive. In view 
of the very interesting consequences which could arise if it turns out to be a first 
order (discontinuous) transition [14], extensive investigations have been carried 
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out into this latter picture. The work reported here is within this scenario in the 
case where the transition starts with the nucleation of hadronic bubbles within 
a slightly supercooled quark-gluon plasma. The bubbles of the new phase then 
grow (the hadron phase is thermodynamically favoured) and subsequently coalesce 
until eventually the universe is filled with the new phase except for remaining dis- 
connected droplets of the quark-gluon plasma. These then shrink and can either 
vanish completely, possibly leaving a significant inhomogeneity in the baryon num- 
ber density [15], or perhaps attain a stable configuration in a new hypothetical 
ground state for the strongly interacting matter {strange quark matter) [16]. 

The present paper is associated with an ongoing research programme 
[9, 10, 11, 17] which has so far focussed on studying the hydrodynamics of the 
growth of single spherical hadronic bubbles within a surrounding quark-gluon 
plasma. While the phase transition directly involves only the strongly interacting 
particles, a crucial role is also played by electromagnetically interacting particles 
(photons, electrons, muons and their antiparticles) and particles which interact 
only weakly (neutrinos and antineutrinos) . These can couple to the quark and 
hadronic plasmas and provide a mechanism for long-range transport of energy 
and momentum on account of having longer mean-free-paths than the strongly 
interacting particles. The effect of this starts to be significant when the radius of 
the growing bubble becomes comparable with the relevant mean-free-path (~ 10 3 
fermi for the electromagnetic interaction and ~ 1 cm for the weak interaction) and 
eventually produces an effective total coupling between the respective components 
when the radius of the bubble is much greater than the relevant mean-free-path. 

Because of the large difference between these mean-free-paths, long-range 
transport can become important at two different stages during the growth of a 
bubble. However, the effects on the hydrodynamics will be exactly the same 
in each case and so we will omit any distinction between them here. The term 
"standard fluid" will always be taken to refer to all of those particles having mean- 
free-paths small compared with the current relevant scale-length for changes in the 
flow of the strongly-interacting matter. 

As mentioned in Section 4, the hydrodynamical equations need to be supple- 
mented by equations of state for each of the two phases. For simplicity, we have 
taken the equation of state for the hadronic matter to be that for an ideal gas of 
massless point-like pions 





1 



(30) 
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where T is the temperature and g is the degeneracy factor (the subscripts h and q 
indicate quantities in the hadron and quark phases respectively). The deconfined 
quarks and gluons cannot be considered as entirely free and for these we have used 
the phenomenological expression given by the M.I.T. Bag Model [18] 



where a positive constant (the "Bag" constant B) is added to the energy density 
and subtracted from the pressure so as to take into account the complex effects of 
confinement. When the photons and relativistic leptons are completely coupled to 
the strongly interacting matter, their contribution can be included by incrementing 



A final comment should be made about the expression adopted in our cal- 
culations for the term (sq) sc , which, as mentioned in Section 3, represents the 
transfer of energy by means of scatterings. While complete expressions for this 
are available for simpler applications [4], the lack of detailed knowledge of the 
interaction processes in the present context has led us to use a phenomenological 
approach, expressing (so) sc as equal to the absorption and emission term, (i.e. 
(e — u>o)/X), multiplied by an adjustable coefficient ranging between zero and one. 
Fortunately, the results of our calculations (which will be presented in paper [11]) 
turn out not to depend sensitively on the value chosen for this coefficient. 

6. Treatment of the Phase Interface 

As long as the radius of the bubble is large compared with the strong- 
interaction length scale, it is reasonable to treat the interface between the hadron 
and quark phases as an exact discontinuity surface with the variables on either side 
of it being linked by junction conditions. The phase interface could, in principle, 
move either supersonically (as a detonation front) or subsonically (as a deflagra- 
tion front) with respect to the medium ahead. In practice, however, the bubbles 
are almost certain to expand subsonically in the present situation (see [9, 19]) and 
we restrict our attention to this case here. 

In contrast with the situation for calculations of detonation fronts, for a de- 
flagration it is necessary to supplement the hydrodynamical equations with an in- 
dependent expression for the rate at which quark matter is converted into hadronic 





(31) 



g h and g q by the relevant number of additional degrees of freedom. 
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matter, derived from considerations of the physical processes occurring at the in- 
terface. A simple expression for this rate was presented in [9] and setting this equal 
to the hydrodynamical flux across the interface, gives the additional equation 



where p s is the interface location, fi s = d\x s jdt and a is an accommodation 
coefficient (0 < a < 1). 

Also, it is very important to pay attention to a correct calculation of the 
causal structure in the vicinity of the interface and it seems that the only satisfac- 
tory way of accomplishing this, when the interface is treated as a discontinuity, is 
by making use of a characteristic method [20]. This involves rewriting the system 
of partial differential equations as a system of ordinary differential equations along 
characteristic curves in the space-time (which can be physically interpreted as the 
worldlines of sonic perturbations in the standard fluid and radiation fluid and the 
flowlines of the standard fluid). The computational technique employed is then to 
use the continuum equations of the previous sections for the bulk of each phase 
and to track the interface continuously through the finite difference grid using a 
characteristic method together with junction conditions (see [10, 11]). 

For each two moment equations which are retained from the infinite hierar- 
chy, there are two families of corresponding characteristic curves with associated 
characteristic speeds. While including a larger number of moments in general in- 
creases accuracy, the role and relevance of the speeds associated with moments 
beyond the first two is controversial. We note that since we are using only the first 
two moment equations in the present work, these difficulties do not arise here. 

In order to write the equations in a characteristic form, it is convenient to 
make use of the following equalities coming from the constraint equations obtained 
from (17) and from the conservation of momentum 




(32) 




4nGabR 



(33) 



Wl, 



a 



(34) 



a 



pw 



and then to rewrite equations (15), (16), (21), (23), (24) as 



aT 



+ B = 0, 



(35) 



u, t + 



bpw 



12 



c„apw „ 
p,t + -^-« lM + B 1 = 0, 



P,t + kf*V + s 2 = 0, 



ct a ( 4 \ 2a 

(u>o),t + + ( 3 + f E ) w ou,n ~ + s 3 = 0, 



(36) 
(37) 

(38) 



(wi), t + ^ Q + (w ) )A1 - ^ Q + f^jwop^ + ^ w i u ,n + 5 4 = 0, (39) 



where 



S = a^GR 

Bi = ac 2 „ 



P + ( 3 + f E ) w o 



gm r 

it 2 



(2u 4ttGR 



-Wi 



,'2u 4ttGR 
= aP ' ~R 



B3 = 2a ^R-a^ Sl ) Wl+ R{3-^) W0 f 



(40) 
(41) 
(42) 



a / 4 



#4 = — — I ^ + / B wosx + 2o 
piy \3 



u _ AnGR 
_ 



r ^iJwi + ^/ B «;o-asi + ^o(/ E ),M- 

(44) 

In deriving equations (35) - (39), the relations = au and _R M = bT have 
been used and we have introduced the local sound speed in the standard fluid 
c s = (dp/de) 1 / 2 . If we now define the state vector 

( U \ 

V 

U= p , (45) 
w 
\wi/ 

equations (35) - (39) can be written in the symbolic form 

£ + < + B-0, ,40) 

where B is the vector whose components are given by (40) - (44) . 

13 



If the expression chosen for the Eddington factor f E is dependent on com- 
ponents of the state vector U, it is necessary to rewrite the partial derivative of 
f E in (44) in terms of the derivatives of the component variables. In doing this 
the elements of the matrix A are obviously modified. For our specific application, 
the expression which we have chosen for the Eddington factor is the following (see 
Appendix B for a discussion of this) 



Equation (39) then becomes 



8u 2 /9 



\h, q 



(l+4 W 2 /3) \Xh,q + R/ ' 



(47) 



a fl \ a { A \ 2a 



with 



a (A \ ( u 4wGR 

B, = -— + f E jw 0Sl + 2a - — Wl ] Wl - 



3aF aTf E w 
+ ^^J E w -as 1 - 



R 



^h,q + R 

(49) 



where, for compactness, we have defined 



K = f E 



Tw Q 



u(l + 4u 2 /3)w 1 
The matrix A then takes the form 

/ aT/bpw 

c 2 s apw/bY 
A = ap/bT 

aw (4/3 + f E )/bT -2aw x jbpw 



(50) 










\ 




a/b 



\ 2aw 1 (l + K)/bT -aw (4/3 + f E )/bpw a(l/S + f E )/b 0/ 

(51) 

Next, we introduce U, the set of left eigenvectors of A, and Aj, the corresponding 
eigenvalues satisfying the relations 



liA Aj/j. 
14 



(52) 



Equation (52) has five distinct eigenvalues (the system is hyperbolic) 



A = 0, (53) 
Ai >2 = ± ^c s , (54) 



A 3 ,4 = ±^^ + / E , (55) 



to which correspond the five eigenvectors 



l = k 0, ± 1, 0, 

d w 



/,._^M±^, i, o, o, o). 



(56) 
(57) 



, . 2(cl-l-K)(l/3 + f E )^ , 1/4 

/ 3 , 4 = M ± r(c 2 _ 1/3 _ /b) -i + f U + /s 



2(/g-2/3-jQ . , /i , 

— j-rWl, 0, ± - + / £ 

pw{c* - 1/3 -f E ) \3 



1/2 



(58) 

where is an arbitrary constant. 

Equation (46) can then be multiplied on the left by U so as to obtain the 
symbolic expression 



dt 1 d\i 



+ liB = 0. 



(59) 



in which each component involves derivatives only along the characteristic direc- 
tion given by d/i/dt = \ . Writing out system (59) explicitly we then get the 
following characteristic form of the equations (35) - (39) 



r f r 

du ± dp + a< — (si ± c s so) + 

pwc s l pw 



+ 4ttGR 



( ^ • /, ) «'o =F r,H';i 



GM 2Tuc s , , 
+ ± — ^ \dt = 0, 



R 2 



R 



(60) 
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which are to be solved along the forward and backward characteristics of the 
standard fluid dp: = ± (a/b)c s dt , 



1 



dw 1 ± [ - + f E ) dw + 



1/2 



+ 



+ 



X 



+ 



X 



3 ' JE 



+ f E )w ± 



2(cg-l-X)(l/3 + /jV 2 
c 2 - 1/3 - f E 



-du+ 



w\dp + a< 



2(f E - 2/3 - K) - 
pw(c* -1/3- f E )_ 

2[(l/3 + f E )(c 2 -l)-Kc 2 ] 



2u 4ttGRw 1 \ 

r r 



R 



c 2 8 -i/s-f E 



Wl ± - + f\ - + fi 



1/2 



w 



+ 



- + f E )W ± 



X 



2(c 2 ; -l-if)(l/3 + /jV 2 
c 2 ~ 1/3 - f E 



Wi 



G 



Ku(l + 4u 2 /3) 

^h,q(l + R/^h,q) 



Wi + 



^3/ B 



+ 



+ 



2c 2 



1/2 



u-T 



Wq 



±1 3+/- 



1/2 



r-u 



pw(cl - 1/3 - /J 

± ? 

. pw(c* - 1/3 -f E ) 



f °-l- K ) 



wi =F ( 3 + / B 



1/2- 



SO + 



1/2 



si = 0, 



(61) 

which are to be solved along the forward and backward characteristics of the 
radiation fluid dp = ± (a/6)(l/3 + f E ) l / 2 dt , and 



dp 5 — dp -dt = 0, 

w 



c 2 s w 



(62) 



which is an advective equation and is to be solved along the flowlines of the stan- 
dard fluid dp, = 0. Finally, R and M are calculated from advective equations 



dR = audi, 



dM = -4iiR 2 au 



dt, 



(63) 
(64) 



16 



and the metric coefficient a is calculated from (25) which is a constraint equation 
on the constant t hypersurface (i.e. it is to be integrated along the direction 
dt = 0). 

The configuration of characteristic curves adjacent to the interface is shown 
in Figure 1 for evolution of the system from time level t to a subsequent time level 
t + At. The dashed lines represent the forward and backward characteristics for 
the radiation fluid r, the full narrow lines are the equivalent characteristics for 
the standard fluids f, the vertical dotted line is the advective characteristic for 
strongly interacting matter in the quark phase and the heavy line is the worldline 
of the interface. 



270 

Figure 1. The configuration of characteristic curves near the phase interface drawn 
in the Lagrangian coordinate framer 



The difference between the s 
the radiation fluid ((1/3 + /J 1/2 



)und speeds in the standard fluid (c s ) and in 
3 large when the former is non-relativistic,but,' 
in the present case, the standard luid is relativistic (with c & — >-~l~/\/3) and the 

is frequent ly-smalk Thisjgarls-te-soms serious 
uations- which -we will -be discussing 



difference between the sound speed; 
complications in numerical solutio: 
in paper [11]. 

The junction conditions linking the values of quantities on gither sic 
interface need to take account of tl e surface tension and surface energy within itr 
A complete relativistic treatment c f this problem can be given most conveniently 
by using the Gauss-Codazzi formalism [21, 22]. In view of some recent comments 
[23] , it is worth stressing that this is an economical way to proceed when working 
within a relativistic Lagrangian framework even in situations where gravity can 
be neglected. 
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Following the approach described in [9] , we can write the following equations 
expressing continuity across the interface of R, dR/dt and ds 



[Rf = o, 



(65) 



[aw + 6/i s T] ± = 0, 



(66) 

[a 2 -6V s ] ± = 0, (67) 

where [A] ± = A + — A~ , {A} ± = A + + A~ and the superscripts ± indicate 
quantities immediately ahead of and behind the interface. 

The junction conditions for the energy and momentum of the standard fluids 

are 

[(e + p)ab] ± = 0, (68) 



[eh 2 \i\ +pa' 1 



af 2 \ 1 d fb 2 (, s \ f M , 2 

T" ) + ^ + jR {b ^ U + aT) 



(69) 



2 [ah dt 

where a is the surface tension (taken here to be independent of temperature) and 
/ = (a 2 — b 2 [jl 2 s ) 1 / 2 . Since the thickness of the interface is much smaller than the 
mean-free-path of the particles of the radiation fluid, it is reasonable to neglect 
any interaction of the latter with the interface, so that the energy and momentum 
junction conditions for the radiation fluid simply reduce to the continuity equations 



abfi s 



+ f E )w - (a 2 + b 2 fi 2 s ) Wl 



0, 



(70) 



a 



^ + f E ) + Irfr \ «•(> - 2«/;// < «';i 



= 0. 



(71) 



The mass function M receives a contribution from the surface energy. At the time 
of nucleation of the bubble, conditions are essentially Newtonian so that 



[M] ± = AnR 2 a, 



(72) 



and the subsequent time evolution is given by 
d 



dt 



[M\~ 



4tvR z 



bT/l s je + w + f w i} - au |p+ Q + f^j w + 



u 



-Wi 



-I ± 



(73) 

Adding conditions (65) - (73) to equations (14) - (16), (21) - (31), and using 
the characteristic method (equations (60) - (64)) to give a correct treatment of 
the causal structure, a satisfactory calculation of the growth of a cosmological 
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hadronic bubble in the presence of long range energy and momentum transfer can 
then be made. 



7. Numerical Computations of Bubble Growth 

A computer code has been constructed in order to implement the formalism 
presented in the previous sections so as to make calculations of the growth of 
hadronic bubbles during the cosmological Quark-Hadron transition. As mentioned 
previously, the code advances the solution from one time- level to the next using a 
standard Lagrangian finite- difference form of the continuum equations for the bulk 
of each phase, while the interface is tracked continuously through the grid using 
a characteristic method together with imposition of junction conditions. This 
code is a further development of the one described in reference [10]. A complete 
account of the numerical techniques used for these new calculations and of the 
results obtained will be presented in our forthcoming paper [11]. 




shows the change irQhe profile of w gs the bubble expands. The horizontal dashed 
line shows the initial profile at the time of b^bjcjtejnu^a^ipjj^ 

As a sample of the results, Figure 2 shows the progressive coupling together 
of the radiation fluid and the standard fluid as a hadronic bubble expands with the 
mean-free-path of the radiation particles becoming progressively smaller compared 
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with the radius of the bubble. Eventually, the coupling becomes essentially total, 
the radiative transfer of energy and momentum ceases to operate and the flow 
tends towards self-similarity. 

8. Conclusion 

Many astrophysical and cosmological situations involve non-stationary spher- 
ical relativistic flows in which a net transfer of energy and momentum occurs be- 
tween a radiation fluid and a standard fluid. In order to describe this kind of 
flow, we have presented here a new set of Lagrangian hydrodynamical equations 
in which the radiative transfer problem is solved using the PSTF tensor formal- 
ism truncated at the second order. The equations obtained in this way describe 
the behaviour of the radiation variables and need to be coupled with the set of 
hydrodynamical equations derived from standard conservation laws. 

One application of particular interest to us concerns bubble growth during 
the cosmological Quark-Hadron transition. We have used our set of equations in 
a numerical study of this problem, investigating the effects of energy and momen- 
tum transfer by weakly and electromagnetically interacting particles. In doing this 
we have treated the interface between the two phases of the standard fluid as an 
exact discontinuity surface and have tracked it continuously through the finite- 
difference grid using a characteristic method together with imposition of junction 
conditions derived with the aid of the Gauss-Codazzi equations. The interface 
moves subsonically and, within our scheme, use of the characteristic equations is 
then particularly important in order to give a correct treatment of the causal struc- 
ture. A computer code has been constructed for solving this system of equations 
and a full account of the techniques used and results obtained will be presented in 
a subsequent paper [11]. 

The set of equations presented here is general in nature and it can be applied 
directly to other situations of interest. Work is now in progress on applying the 
equations of Sections 3 and 4 to the study of non-stationary spherical accretion 
onto a black hole. A computer code has been constructed for this and is at present 
under test. 
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Appendix A 

In this Appendix we give a brief sketch of the calculations leading from 
equations (12) - (13) and (17) - (20) to the hydrodynamical equations (15), (16) 
and (21) - (29). For clarity, we make a separation into two sub-sections, one for 
the radiation fluid alone and the other for the combined fluids. 



THE RADIATION FLUID 

As mentioned in Section 3, we use the spherically symmetric line element (5) 
and expression (4) for the stress-energy tensor of the radiation fluid. The first three 
PSTF moments are written as 

M = w , (Al) 
M a = w 1 e%, (A2) 

and 

= w 2 (e?e? - - \e%e^ , (A3) 

where wq, w\ and w 2 can be physically interpreted as explained in Section 3 and 
(e 5 , e f , eg, e^) is the orthonormal tetrad carried by an observer comoving with 
the standard fluid. The components of this tetrad are 

= Q, 0, 0, o), (A4) 
e?=(o, \, 0, o), (A5) 

e * = (A6) 

e H°'°'°'^)- (A7) 
The comoving observer's four- velocity u and four-acceleration g are given by 

u a = (-, 0, 0, o] =e?, (A8) 



a 



and 



S a =(o, 0, o), (A9) 
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(so that = g a /g, where g = (g 0l g i) 1 ^ 2 )- The covariant derivative of the radiation 
fluid stress energy tensor is 

Tf = i (m.^u 13 + MvF.puP + Mu a v?.^j + 

+M a . f3 u f3 + M a u p . p + M p . )f3 u a + M p u a . p + M ap . p , 

(A10) 

where © = u a . a is the expansion and u a .pV,P = g a . The contraction of expression 
(A10) with u a then gives 

««T°% = -M,pvP - - 9 — Wl ^+ 
' 6 g 

+wi(— \ vPu a - wi ( — ^ +M af3 . a u a . 

(All) 

After some further manipulation of expression (All) and using the equality 



+W 2 



(A12) 

it is then possible to rewrite equation (12) in the form (15). We proceed in the 
same way with the derivation of equation (13). Bearing in mind that njP^ a = S 1 ^ 
it follows from (A10) that 



„i 

+M\ p u fi + w 1 j<?> + M p u\p + M^. p . 



(A13) 

Writing out each of the terms explicitly and using the expression 

^,=v>^M^ +3 ^y <ai4) 

equation (A13) can finally be recast in the form (16). 
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THE COMBINED FLUIDS 



Expressions for the non-zero components of the total stress-energy tensor 
(T a/3 = T^P + T^P ) can be calculated using equation (4) and expressions (A4) - 
(A7), and this gives 

— (e + w ), 



cr 

ol _ wi 
ab ' 
1 



jiOO 

r 

rpll _ 

rp 22 

rp33 _ 



W 

b2 yP+ T +^ 



1 ( Wq W 2 

R 2 sm 2 9V 3 2 '* 



(A15) 
(A16) 

(A17) 
(A18) 
(A19) 



where it is easy to see the contributions coming from the standard fluid and the 
radiation fluid respectively. 

We will not write down here the expressions for the non-zero Christoffel 
symbols and the relevant components of the Ricci tensor, which are obtained by 
straightforward but particularly tedious calculations, but proceed directly to the 
form of the four independent Einstein equations 



(T° ) 

8nG(e + w )R 2 R yll + 



8nGbR 2 R 



a 



-W! 



R 



a 



R 



+ 1 



(A20) 



W 



8wG\p+-£- + w 2 1 + 



8nGaR 2 R 



-Wi 



-u 



Rt 

a 



R 



+ 1 



(A21) 



( 4 
8txG\v + e + -wq - 2w 2 



(T°) 



8TVG-W! = 

a 



1 

ab 



a 



+ 



1 

+ R 2 



Rt 

a 



i 2 R 



tfl 



R 



+ 1 



(A22) 



(A23) 
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(the symbol in brackets identifies the component referred to). Equation (A23) is a 
constraint equation which, in the form 



b 4 = -4~ ( R ^ - — R J - ^GabR Wl ) 
o it, M V a J 



(A24) 



has been used frequently in the course of the calculations outlined in this Appendix. 
In particular, writing out explicitly the time derivative on the right-hand side of 
equation (A21) and making use of the (A24), the expression 



AnGRawi — 



R,t ( p ,,j, + bsi 
e + p 



(A25) 



can be obtained and this can then be further transformed, by means of equation 
(27), so as to arrive at the form (21). 

Next we turn to writing out explicitly the hydrodynamic conservation 
equations (18) - (20) which then take the form 



-U r 



T a<3 , = = s + (e + p)^vP + (e + p)6 - p,0 



a/3 



(A26) 



njP^T^.p = = - 8l + (e + p)n 3 v?.^ + (g af3 + u a u p )p^ , (A27) 



{pu a ). a = = p yCt u a + pu a a + pu c 



(y/det(-g^)) a 
v /det(-(7«/ 3 ) 



(A28) 



Using the expressions V = R^/b and u = Rj/a together with (A24), equations 
(A26) - (A28) can be converted to the final form given in equations (23) - (25). 
Finally, we note that if we rewrite equations (A20) and (A21) as 

{ 2! [(^) 2 " (^) 2 + 1 }}^= ^ R2R ^ + + = M,„ (A29) 



and 



R 
2G 



a 



R 



+ 1 



,t 



-InR- Rj ( p+ -w + +w 2 ) = M,t- 



(A30) 

this gives expressions (26) and (29) for the generalized mass function M. These 
reduce to the familiar expressions for the standard mass function when the 
radiation terms are omitted. 
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Appendix B 



As mentioned in Section 3, a particularly important point in the PSTF 
approach is the truncation of the infinite hierarchy of moment equations by means 
of a suitably defined closure relation derived on the basis of physical considerations. 
We have used the closure equation (14) and, for our particular problem, an 
appropriate expression for f E is given by equation (47). This appendix gives our 
justification for making this choice. 

A key point in choosing an expression for f E is that it should have the 
correct asymptotic behaviour in the optically-thin and optically-thick limits. We 
are considering the case of a single spherical hadronic bubble which is initially 
nucleated at rest with a radius small compared with the mean-free-path of the 
radiation. Under these circumstances, the radiation field will be everywhere rather 
accurately uniform and isotropic (unless there is some other perturbing influence) 
and since the bubble radius is very small compared with the horizon scale, it is 
also a good approximation to neglect cosmological expansion. Since wo, w\ and 
«)2 are all measured with respect to the local rest-frames of the standard fluid, the 
values which they take during the early part of the bubble expansion are those 
produced by motion of the fluid rest frames with respect to the essentially uniform 
radiation field. 

When the radiation is isotropic in its mean rest frame, its stress-energy tensor 
takes the perfect fluid form 



where e R and p R are the radiation energy density and pressure (p R = e K /3), 
measured in the mean rest frame of the radiation, and v a = dx a jdr is the four- 
velocity of this frame relative to some specific observer. For purely radial motion 
in our metric 



Tf = (e R +P R )v a vP+p R g 



(Bl) 





and the non-zero components of four- velocity are 




(B3) 



where 



7 = 



l 



and 




(B4) 
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To find the value of fi for the radiation frame with respect to the standard fluid, 
we note that since each element of the radiation fluid is remaining at a constant 
value of R 

dR = R it dt + R^dfi = audi + bTdfi = 0, (B5) 



and so 



^ = ~5f' (B6) 



which leads to the following expressions 



u 



« = - F ^ = rw- (B7) 

If we now write the stress-energy tensor (Bl) in the frame comoving with the 
standard fluid, we can compare the new expressions for the components with the 
ones appearing in equation (4) and so obtain the following system of equations 

i)+r*(-h)' (B8) 
V) • < B9 ) 



rpOO _ 


W 


4 


R 


a 2 


" 3 


rpOl 


Wi 


4 


R 


ab 


" 3 


Jill 


&2 ( 


'wo 


R 




v 3 



a 

2 



r-A , 1 1 



- 3 e » ) + 3 C » U ) ■ (B10> 

The solution of this system then leads to the expressions 

w = ^-(3 + v 2 )e R , (Bll) 

= ^7 2 ^e H , (B12) 

w 2 = £(l + ti*)e R -^. (B13) 

If we now define (wo) N = e R to be the radiation energy density at the bubble 
nucleation time (when there is no fluid motion), equations (Bll) - (B13) can be 
suitably transformed so as to give expressions for the energy density, flux and shear 
scalar of the radiation as seen from the standard fluid 

wo=(l + \ V 2 U _ U ? ) (™oL - (l + \u 2 ^j (w ) N , (B14) 

4 Tu 4 
Wl = ~3 T 2 -u 2 ( Wo ^ N ~ ~3 Tu ( W( >)»> ^ B15 ^ 

8 u 2 8 
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The approximate forms of the expressions (B14) - (B16) result from noting that 
since the dimensions of the bubble are small compared with the horizon scale, 

r 2 - u 2 ~ 1. 

Equations (B14) and (B16) together with the definition (14), give the 
following analytic expression for the variable Eddington factor during the first 
stages of the bubble expansion 

f 8u2/9 mi7i 

This is the "optically thin" limit. At the other extreme, the "optically thick" limit 
arises when the radius of the bubble is large compared with the radiation mean- 
free-path and complete coupling has been attained between the radiation and the 
standard fluid over length-scales comparable with the radius of the bubble. When 
this happens, interactions make the radiation isotropic in the local fluid rest frame 
so that W2 — > and f E — > 0. A suitable smooth join is required in between the 
two asymptotic limits and to do this we have multiplied the expression in (B17) 
by Xh, q /(^h,q + R) which then gives equation (47). Experiment has shown that 
reasonable variation in the form of the join makes an insignificant change in the 
results obtained. 
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